Scale-free networks are not robust under neutral evolution 
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Recently it has been shown that a large variety of different networks have power-law (scale-free) 
distributions of connectivities. We investigate the robustness of such a distribution in discrete 
threshold networks under neutral evolution. The guiding principle for this is robustness in the 
resulting phenotype. The numerical results show that a power-law distribution is not stable under 
such an evolution, and the network approaches a homogeneous form where the overall distribution 
of connectivities is given by a Poisson distribution. 
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There are many different areas of contemporary science 
where the concept of networks are of special importance. 
Recently, the amazing result that such diverse networks 
as the World Wide Web |Q , collaborations of movie ac- 
tors Q , the electrical power grid of western USA || , cita- 
tion patterns of scientific publications Q and metabolic 
networks j| , all have distributions of connectivities that 
are scale- free, i.e., obey some power- law. One reason for 
this quite general behaviour seems to be that the net- 
works are grown by addition of new nodes and that each 
such new node preferentially attaches to other nodes with 
high number of connections These networks are also 
more robust against unintelligent attacks than homoge- 
neous networks, where each node has approximately the 
same number of connections ||. For the case of genetic 
regulatory networks, mutation is an example of such an 
"unintelligent attack", and it is perhaps not surprising 
that the same type of scale-free distributions are found 
in the metabolic networks of so far 43 different organisms 
from all three domains of life (bacteria, eukarya, and ar- 
chaea) B. 

The evolution of life is a random process with selection 
H , although all details about how this occur with inter- 
actions among genotypes, phenotypes, and environment 
are not totally clear. Neutral evolution is the hypothesis 
that evolution mainly proceeds as a random walk which 
does not affect the phenotypes Qj. Experimentally, it is 
supported on the microlevel by the fact that most of the 
important macromolecules of life have forms which are 
functionally identical variants. If this idea of neutrality 
is correct also on a higher level, it means that the useful- 
ness of fitness landscapes for describing the evolution of 
life as a hill climbing process is limited. 

Here we explore the idea of neutral evolution and how 
the distribution of connectivities in a genetic regulatory 
network changes under mutations that are phenotypically 
silent. The fundamental constituents in our model are 
the genes of the organism, represented by the nodes of 
the network. It has been suggested that such a system 



can be well approximated by a Boolean network ||, be- 
cause of the "on-off" nature of the biochemical switches. 
The exploration is performed by simulating evolution in 
discrete threshold networks with robustness as the guid- 
ing principle for when a mutation will survive to future 
generations. We find by numerical simulations that the 
scale-free distribution cannot be maintained under neu- 
tral evolution in such networks. Instead, the networks 
evolve towards a Poisson distribution, regardless of the 
actual realization of the initial scale-free distribution. 

The discrete threshold network is composed of TV 
nodes, Oj, which are connected by a square matrix with 
elements Wij. The values of the nodes are cr, G { — 1,1}, 
representing the corresponding gene to be expressed (+1) 
or not (—1). The coupling matrix takes values £ 
{—1,0,1}, with +1 if gene j is an activator of gene i, 
— 1 if it is a repressor, and if no connection exist. The 
dynamics of the network is described by the updating 
rule 



Oi(t + 1) = sgn 
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where the sign-function is defined as —1 for all negative 
arguments, and +1 otherwise (including zero). Of special 
importance is the mean number of connections to each 
node, K, which is calculated as 
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that is, we make no distinction between repressors and 
activators, and it has the same value for both ingoing 
and outgoing connections. 

This is a special case of Boolean networks, with similar 
structural and statistical properties || . These properties 
include both transients and limit cycles (attractors), as 
well as phase transitions for a specific critical connec- 
tivity, K = K c . An analytical approach is limited due 
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to the non-Hamiltonian character of the system, but re- 
sults within the so-called annealed approximation show, 
for Boolean networks, that for K below the critical con- 
nectivity K c = 2, there are many disconnected regions, 
while above K c most of the nodes are connected and the 
limit cycle period increases exponentially with the num- 
ber of nodes. Also, at least in some intervals above K c , 
the size of the attractors diverges almost exponentially 
with increasing connectivity |L0]]. Note here that we use 
a form where the number of ingoing connections might 
differ from the number of outgoing, i.e., the matrix does 
not have to be symmetric, and that we do not impose 
the restriction that the number of connections should be 
the same for all nodes ||. 

If every two nodes in a network are connected with the 
same probability p, we have the Erdos-Renyi model for 
random graphs [p"l| , sometimes referred to as a homoge- 
neous network. It is well known that in such networks, 
the number of connections to each node follows a Poisson 
distribution (with exponential decaying tails), and hence 
there will hardly be any nodes with a large number of 
connections. Loosely speaking, this is the most common 
form of a random graph. 

To simulate neutral evolution, we start by generating 
a network with elements Wij by the procedure described 
in B, i.e., we start from a small random network, and 
add new nodes by preferential attachment. The result is 
a scale-free network with a probability for a given node 
to have K links (either ingoing or outgoing) proportional 
to K^ 1 for K > 2, with, in our case, 7 « 1.63 (solid line 
in Fig. |l|). We use consequently in this paper N = 1024, 
which results in an average connectivity value initially 
slightly below the critical value K c = 2. Interestingly, 
a recent letter showed that another form for evolving a 
discrete threshold network (adding links to quiet nodes, 
removing links from active) leads to an average connec- 
tivity of 2.55 for this size of network [O. The sign of a 
specific connection specifies if we have an activator (pos- 
itive value) or a repressor (negative value). These signs 
are here chosen randomly with equal probabilities. 

The evolution now proceeds by the following proce- 
dure: The network is mutated by either 

1. One non-zero element is put to zero (connection 
removed) 

2. One zero element is turned into ±1 (with equal 
probability) (connection added) 

3. Both of the above, i.e., one connection is added and 
another removed 

These three alternatives (in the order given) occur with 
the probabilities 0.300 : 0.333 : 0.367, forming a new, 
mutated network. The values of these probabilities were 
chosen to obtain a network with a relatively constant 
mean number of connections also for the comparatively 



small number of generations and the initially low number 
of mean connectivity we consider. However, also other 
values have been tried, and the results do not depend 
critically on their exact magnitudes. To either reject or 
accept the new, mutated network, we use robustness as 
the guiding principle. This is achieved by picking by ran- 
dom an initial state, {ci}, with equal probability for each 
single node Ci being either positive or negative. This 
state is iterated in both the original and the mutated 
network until we either enter into the same limit cycle 
in both networks, or the two iterants cease to coincide. 
In the former case we accept the mutated network, and 
replace the original one with the mutated version. This 
is then an evolutionary step within the neutral evolution. 
In the latter case, we reject the new, mutated network, 
since the effect of the mutation were not silent. Finally, 
we return to the mutation step and repeat the procedure. 
Notice this introduces two different time scales in the evo- 
lution. The one corresponding to the iteration of states 
{o~i} relates to a single generation, while the much slower 
process of accepting or rejecting new networks, i.e., the 
rewiring of connections toy , corresponds to the evolution 
over generations. 

This way of simulating neutral evolution has earlier 
been explored by Bornholdt and Sneppen, both by truly 
Boolean networks |l^] and by discrete threshold networks 
]l4| |. They studied, however, the phenomenon of punctu- 
ated equilibrium and distribution of waiting times, and 
ignored the distribution of connections. Their study clar- 
ified that this model exhibit many of the known prop- 
erties of evolution, such as 1// power spectra |l5| and 
1/t 2 stability distribution |l6| ], in accordance with simi- 
lar scalings found in the statistics of birth and death in 
the evolutionary record. 

In Fig. ^ we show the initial distributions of con- 
nections and the distribution after 30 000 generations 
in one evolutionary run for the number of connections 
leading into the nodes. The result for connections lead- 
ing out from them are quite similar, and are for clarity 
not shown. Although the limited number of nodes (due 
to computational constraints) makes the statistics some- 
what fuzzy, it is still clear that the distribution changes 
from a power-law to an approximately Poisson distribu- 
tion. To get a better picture, we have used the well- 
established technique of binning the values for the initial 
distribution. The solid and dotted lines are a power-law 
and a Poisson distribution, respectively. The exponent of 
the straight line P{K) ~ K 1 is found by a least squares 
fit to be 7 ks 1.63. This should be compared to the the- 
oretical value given in Q for a fully directed graph of 
infinite size, which is 7 = 2. Although the number of 
nodes we use is small compared to the networks consid- 
ered there, the correspondence seems acceptable. The 
Poisson distribution drawn is for the expectation value 
estimated by the mean value of connectivities, K, at the 
actual generation. This means that in the general Pois- 
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son distribution function 

P(iO = ^exp(- M ), (3) 

we estimate the parameter \i with the mean connectivity 
K. No curve fitting is used this time, but nevertheless 
the correspondence is quite striking. 

In Fig. ||, we show the variation of mean connectivity, 
K, for the first 30 000 generations. Because of the defini- 
tion (0), this is the same value both for outgoing and in- 
going connections. The curve shows that the mean num- 
ber of connections remains fairly constant, although the 
detailed dynamics is non-trivial (with punctuated equi- 
libria, etc., as discussed in (l^jlj]). It also shows that 
we in this run constantly are below the critical value of 
K c = 2. This is, however, not a critical aspect of the sim- 
ulation, which other runs (not shown) have indicated. 
Hence any change in distribution among these connec- 
tions cannot be due to changes in the mean connectivity. 
Nevertheless, there is according to Fig. [j] a real change 
in distribution from the start of the simulated evolution 
to the end of our calculations. To shed some further 
light onto the transition from a scale-free network with 
a power-law distribution to a homogeneous network with 
a Poisson distribution, we calculate the weighted mean 
square deviation 

1 N 

d 2 = -J2 K H K )- NP ( K )\ 2 > ( 4 ) 

where n(K) is the number of nodes with K connections 
and P{K) is the Poisson distribution (J|). For each com- 
parison, we use the actual value of K as estimate for the 
expectation value /i, but no other fitting is performed. 
Because the tail of the distribution function is the most 
critical, we give higher weight to larger number of connec- 
tions by multiplying each term with the actual number of 
connections. The results are shown in Fig. ||. Although 
the exact details for the ingoing and outgoing connec- 
tions differ, and it is clearly seen that the distributions 
eventually approach the Poisson distribution, regardless 
of which measure we consider. 

The lengths of the limit cycles for the accepted net- 
works in this evolutionary scenario vary between 1 and 
20, with an average of approximately 6. The transients 
have lengths between 10 and 35 steps. Both these results 
indicate that we are in the regime of many, small, dis- 
connected attractors, which is fully consistent with the 
mean connectivity K being less than the critical value 
Kc = 2. 

Better statistics, i.e., less fuzzy distributions, are ob- 
tained if we change the mutation rule above somewhat. 
Instead of having the possibility to separately add or re- 
move a connection, we stick solely to alternative three, 
which means that the number of connections, and hence 



K, remains constant. This is clearly a less realistic sce- 
nario than before from a biological point of view, but can 
help to see what distributions we really have. In Fig. ||, 
we show the mean values of the distributions for all gen- 
erations between number 200 000 and 300 000, when we 
start with the same initial scale-free network used to ob- 
tain the results of Fig. [I]. The dotted line is the theoreti- 
cal Poisson distribution for the actual mean connectivity 
(K = 1.9678). The weighted mean square deviations 
from this theoretical value are 1.18 for the ingoing con- 
nections and 0.06 for the outgoing, respectively. To the 
prize of having incorporated an unrealistic restriction, we 
have obtained distributions which are considerably closer 
to the theory. 

To check the robustness of these results, we have re- 
peated the calculations many times with different forms 
of initial network in the construction of the original scale- 
free network, as well as checked many different realiza- 
tions. We have also started with networks with a power- 
law distributed number of connections obtained directly 
from a random number generator, i.e., without the pro- 
cess described in 0. In all these cases, our results do not 
change, i.e., the systems always end up with a Poisson 
distribution consistent with the random graph theory. 

In conclusion, we have studied the evolution of ini- 
tially scale-free networks, i.e., networks where the distri- 
bution of the number of connections to each node fol- 
lows a power-law. The networks are evolved under the 
hypothesis of neutral evolution, which is implemented 
as a robustness criterium for the limit cycles in discrete 
threshold networks. The result is that the scale- free dis- 
tribution is not robust under such an evolution, but in- 
stead all networks end up in a homogeneous form, with a 
Poisson distribution of connectivities. This result is sur- 
prising, since it has been shown that a scale-free network 
is more robust towards random attacks than an expo- 
nential network § . Although one should be very careful 
with infering any definite statements with respect to bi- 
ology from such simple models as the one presented here, 
we can speculate and draw the tentative conclusion that 
the addition of new nodes with preferential attachment 
seems to be a force that manages to repress the changes 
due to neutral evolution. Future studies might shed some 
light on the presumably different time-scales that are ac- 
tive here. 
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FIG. 1. Distribution functions for number of ingoing con- 
nections for an initially scale-free network with N = 1024 
nodes at two different generations. Diamonds: Initial 
power-law distribution. Crosses: Distribution after 30 000 
generations of neutral evolution. The full line is a power-law, 
P(K) ~ K~"< with 7 = 1.63, and the dotted line is a Poisson 
distribution (see text for details). 



FIG. 3. Weighted mean square deviations from the Poisson 
distribution for an initially scale- free network, evolving under 
neutral evolution, where the expectation value at each gener- 
ation is estimated by the actual mean connectivity. Full line 
represents ingoing connections, dotted line represents outgo- 
ing connections. 
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FIG. 2. Mean connectivity for an initially scale-free net- 
work evolving under neutral evolution. 



FIG. 4. Mean distribution for all generations between 200 
000 and 300 000 for an initially scale-free network evolving un- 
der neutral evolution with constant connectivity. Diamonds 
are ingoing connections and crosses outgoing. The dotted line 
is a Poisson distribution with expectation value estimated as 
the mean connectivity. 



5 



